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Abstract 

We investigate in detail two models describing how stresses propagate and fluctuate 
in granular media. The first one is a scalar model where only the vertical component 
of the stress tensor is considered. In the continuum limit, this model is equivalent to a 
diffusion equation (where the role of time is played by the vertical coordinate) plus a 
randomly varying convection term. We calculate the response and correlation function of 
this model, and discuss several properties, in particular related to the stress distribution 
function. We then turn to the tensorial model, where the basic starting point is a wave 
equation which, in the absence of disorder, leads to a ray-like propagation of stress. In the 
presence of disorder, the rays acquire a diffusive width and the angle of propagation is 
shifted. A striking feature is that the response function becomes negative, which suggests 
that the contact network is mechanically unstable to very weak perturbations. The stress 
correlation function reveals characteristic features related to the ray-like propagation, 
which are absent in the scalar description. Our analytical calculations arc confirmed and 
extended by a numerical analysis of the stochastic wave equation. 



1 Introduction 



Granular media are materials where stress fluctuations are large, even on scales much larger 
than the grain size. Repeatedly pouring the very same amount of powder in a silo results in 
fluctuations of the weight supported by the bottom plate of 20% or more ||. This weight 
furthermore changes very abruptly when temperature changes by only a few °C, which induces 
only very small changes of the size of each grain ^ . 

More quantitative experiments were recently performed by Liu et al. Brockbank et al. 
m and Mueth et al. @, where the local fluctuations of the normal stress deep inside a silo or 
at the base of a sandpile were measured (see also and for early qualitative experiments 
Q). It was found that the stress probability distribution is rather broad (i.e. the relative 
fluctuations are of order one), decaying exponentially for large stresses. A simple 'scalar' 
model for stress propagation was introduced and studied in detail ^, which predicts a 
stress probability distribution in good agreement with experimental (and numerical) data. 
However, this model only considers the vertical normal component of the stress tensor, and 
discards all the other components: in this sense the model is scalar. 

A fully 'tensorial' model for stress propagation in homogeneous granular media was pro- 



posed in [10, 11, |l^ to account for the pressure 'dip' which is observed experimentally below 
the apex of conical sandpiles. The most striking feature of this model is that the stress 
propagation equation is a wave equation, with the vertical axis playing the role of time. Cor- 



respondingly, the stress propagates (in two dimensions, see [|10|) along two rays which makes 



a certain angle with the vertical axis (the 'light cone'). This must be contrasted with the 
scalar model, where stresses travel essentially vertically, which predicts a central pressure 
'hump' (rather than a 'dip'). 

It is thus a priori not obvious that the scalar model is a suitable starting point for the 
description of fluctuations. Conversely, the influence of local randomness within the tensorial 
model was not yet investigated, and is very interesting per se. In particular, it is important 
to know if and how the idea of a 'light cone' survives in the presence of disorder, and how 
the stress fluctuations develop. 

The aim of the present paper is to calculate analytically (in two dimensions) the average 
response function (Green function) and the two-point correlation function for the tensorial 
model in the presence of disorder, and to compare the results to those obtained within a scalar 
description. We find that the cone survives at small disorder (although the cone angle is shifted 
and acquires a non zero width, which we compute). More surprisingly, the Green function 
takes negative values |^, a feature which we checked numerically, and which we discuss in 
detail in terms of the essential "fragility" of the contact network. We show that the two-point 
correlation function keeps a signature of this cone like propagation. For large disorder however, 
the theory suggests that the structure of the large scale equations could change drastically, 
from an hyperbolic wave equation to an elliptic equation, akin to (but distinct from) those 
appearing in elasticity theory. The interpretation of the equations however suggests that by 
the time this happens, the pile is unstable to any perturbations and spontaneously rearranges. 

The 'tensorial' stress probability distribution is investigated numerically, with certain 
results which are close to those of the scalar model. We explain this by showing that a special 
case of the tensorial model actually reduces to the superposition of two independent scalar 
models. 

This paper is constructed as follows: in section 2, we review the properties of the scalar 
model, including results which appeared in the literature in very different contexts (scalar 
diffusion in turbulence, localisation). In section 3, the 'random wave equation' for the tensorial 

^That the Green function can take negative values in the presence of inhomogeneities was aheady noticed 
within the fpa model in . 



case is motivated by a microscopic model and simulations, and studied using perturbation 
theory in the strength of the disorder. We discuss how the line shape of the response function 
distorts from two delta peaks to (eventually) one broad peak as disorder increases. Some 
generalisations of the the 'random wave equation' are considered in section 4. In section 5, 
we present numerical results for the stress distribution function and compare them with the 



predictions of the scalar model, and also of direct simulations of sphere packings [14, 15 1. We 
discuss a limit where the two models can be quantitatively compared. Finally, in section 6, 
a summary of the most interesting results is given, with suggestions of new experiments and 
open questions. 



2 The Scalar Model 
2.1 The discrete version 

• Definition. 

The main assumption of the scalar model is that only the vertical normal component of 
the stress tensor w = Uzz (the 'weight') needs to be considered. If the grains reside on the 




Figure 1: The 'Chicago' model with two neighbours, gi's are independent random variables, except 
for the weight conservation constraint: q+{i,j) + q^{i,j) — 1. 



nodes of a two-dimensional lattice (see figure ||), the simplest model for weight propagation 
down the pile is: 

w{i,j) = Wg + q+{i - 1, j - l)w(i - 1) +g_(i + l,j - l)w{i + - 1) (1) 

where 'wg' is the weight of each grain, and q±{i,j) are 'transmission' coefficients giving 
the fraction of weight which the grain (i,j) transmits to its right (resp. left) neighbour 
immediately below. Mass conservation imposes that g+(z,j) +q-{i,j) = 1 for all i,j's. The 
case of an ordered pile of identical grains would correspond to q± = The authors of [Q, 
proposed to take into account (in a phenomenological way) the local disorder in packing, 
grain sizes and shapes, etc., by choosing q+{i,j) to be independent random numbers (except 
for the above constraint), for example uniformly distributed between and 1. This model, 
which we shall call the 'Chicago' model (or q model), was originally written with an arbitrary 



number N of downward neighbours {N = 2 in the example above), and can thus be (in 
principle) generalized to three dimensions. 



• Results on the stress distribution. Universality ? 



The case of a uniform distribution of the g's is interesting because it leads to an exact 
solution for the local weight distribution P{w). In this limit, the correlation between two 
neighbouring sites at the same altitude j is zero for all j. For more general q distributions, 
this is true only when j is large (see below). Thus P^w) obeys the following mean- field 
equation: 



y[w) = dqidq2p{qi)p{q2) dwidw2Pj{wi)Pj{w2)S[w - {wiqi + W2q2 + wq)] (2) 
Jo Jo 



where p{q) is the distribution of q, here taken to be p{q) = 1. In the limit j oo, the 
stationary distribution P* of this equation is given by: 

^'*H = ^exp-^ (3) 

where 2W = jwg is the average weight. For N ^ 2, the distribution is instead a F distribution 
of parameter N; its small w behaviour is w^~^ while the large w tail is exponential. Liu et 
al. ^ have argued that this behaviour is generic: for example, the condition for the local 
weight w to be small is that all the N q^s reaching this site are themselves small; the phase 
space volume for this is proportional to w^~^ if the distribution p{q) is regular around q = 0. 
However, if instead p{q) oc q'^~^ when q is small, one expects P*{w) to behave for small w 
as w~°', with a = 1 — N'j < 0. Similarly, the exponential tail at large w is sensitive to the 
behaviour of p{q) around g = 1. In particular, if the maximum value of q is qM < 1, one can 
easily show by taking the Laplace transform of equation (|2|) that P* (w) decays faster that 
an exponential: 

a log N 

log P*(w) oc -w^ with /?=— I- (4) 

(Notice that /3 = 1 whenever qM = 1, and that /3 ^ oo when qM = 1/-^)- 

In this sense, the exponential tail of P*{w) is not universal: it requires the possibility that 
one of the q can be arbitrarily close to 1. This implies that all other g's originating from 
that point are close to zero, i.e. that there is a nonzero probability density that one grain is 
entirely bearing on one of its downward neighbours. 

Note that if q can only take the values or 1, the distribution P{w) becomes a power law, 
P*{w) oc 'w~°', with a = 4/3 for N = 2 |^]. This power law is however truncated for large w 
as soon as values for q different from and 1 are allowed. 

How well does the simple distribution (^ compare to experiments and numerical simu- 
lations? The exponential decay for large w appears in some cases to overestimate both the 
experimental and numerical tail |15| (see also section ||), suggesting a value of P somewhat 



larger than 1. On the other hand, the probability to observe very small w is much underesti- 
mated by equation (^): see [^, 14, |^ and section ^. This might be due to the fact that arching 



effects are absent in this scalar model. A generalisation of the Chicago model allowing for 
arching was suggested in , which generates sites where g+ = 1 and g_ = (or vice versa) . 
This indeed leads to much higher probability density for small weights, P*{w) oc w~" as 



argued in - see also | It ] 



2.2 Continuous limit of the scalar model. 



Let us focus on the case N = 2 and define v to be such that q±{i,j) = (1 =b v{i,j))/2. If v is 
small, the local weight is smoothly varying, and the discrete equation (||) can then be written 
in the following differential form: 

dtw + d^{vw) = p + Dod^^w (5) 

where x = ia and t = jr are the horizontal and (downwards) vertical variables corresponding 
to indices i and j of figure |l], and a and r are of the order of the size of the grains. The 
vertical coordinate has been called t for its obvious analogy with time in a diffusion problem. 
p is the density of the material (the gravity g is taken to be equal to 1), and Dq a 'diffusion' 
constant, which depends on the geometry of the lattice on which the discrete model has been 
defined. For a rectangular lattice as shown in figure ||, Dq = More generally, the diffusion 
constant is of the order of magnitude of the size of the grains, a. 

In this model and in the following, we shall assume that the density p is not fiuctuating. 
Density fiuctuations could be easily included; it is however easy to understand that the 
resulting relative fluctuations of the weight at the bottom of the pile decrease with the height 
of the pile H as H~^/'^, and are thus much smaller than those induced by the randomly 
fluctuating direction of propagation, encoded by q (or v), which remain of order 1 as H ^ oo. 

Two interesting quantities to compute are the average 'response' G(x,i|xo,to) to a small 
density change at point (xo,to)) measured at point (x,i), and the correlation function of the 
force field, C{x,t, x' ,t') = {w{x,t)w{x' ,t')) ^ (connected part), where the averaging is taken 
over the realisation of the noise v{x,t). 

Equation (|5[) shows that the scalar model of stress propagation is identical to that de- 
scribing tracer diffusion in a (time dependent) flow v{x, t). This problem has been the subject 
of many recent works in the context of turbulence we believe that interesting qual- 

itative analogies with that field can be made. In particular, 'intermittent' bunching of the 
tracer field correspond in the present context to patches of large stresses, which may induce 
anomalous scaling for higher moments of the stress field correlation function. We refer the 
reader to [17, IS] for further details. 



• Statistics of the noise v{x,t). 



The noise term v represents the effect of local heterogeneities in the granular packing. 
Its mean value is taken to be zero, and its correlation function is chosen for simplicity to 
be of the factorable form {v{x,t)v{x' ,t')) = a'^gx{x — x')gt{t — t'), where gx and gt are noise 
correlation functions along x and t axis. We shall take gx and gt to be short-ranged (although 
this may not be justified: fluctuations in the microstructure of granular media may turn out 
to be long-ranged due to e.g. the presence of long stress paths or arches), with correlation 
lengths £x and if Our aim is to describe the system at a scale L much larger than both the 
lattice and the correlation lengths: ti, T, £x^ i-t This will allow us to look for solutions 

in the regime k,E ^ 0, where k and E are the conjugate variables for x and t respectively, 
in Fourier-Laplace space. However, we shall see below that the limit a,T,ix,it can be 
tricky, and must be treated with care: this is because the noise appears in a multiplicative 
manner in equation (^) 0. For computational purposes, we shall often implicitly assume that 
the probability distribution of v is gaussian; this might however introduce artefacts which we 
discuss. 



■^In the tensorial case, the hmit it ^ actually makes the problem trivial, for a reason which will become 
clear below. 



Fourier transforms. 



The limit where a^lx — > is iU defined and leads to a divergence of the perturbation 
theory in a for large wavevectors k. We thus choose to regularize the problem by working 
within the first Brillouin zone, i.e., we keep all wave vector components within the interval 
X = [—A, +A], where A = Our Fourier conventions for a given quantity / will then be the 
following: 

/•A fii, 

f{x,t) = / -e"^^f{k,t) (6) 
J-A 27r 

f{k,t) = 4 E (7) 

x=—oo 

One has to be particularly careful when computing convolution integrals, such as / ^fi{q)f2ik- 
q) which must be understood with limits —A + /c, A (resp. —A, A + fc) if /c > (resp. k < 0). 
An important example, which will appear in the response function calculations, is: 

% = l!^ + oie, (8) 



L 



Let us then take the Fourier transform of equation (|5|) along x, to obtain: 

{dt + DQk^)wk = pk + ik J ^WgVk-g (9) 

Our aim is to calculate, in the small-A: limit, the average response (or Green) function G{k,t — 
t') defined as the expectation value of the functional derivative {5w{k,t)/5p{k,t')); and the 
two points correlation function of w, {w{k,t)w{k' ,t)) = 27rd{k + k')C{k,t). 



• The noiseless Green function. 



The noiseless (bare) Green function (or 'propagator') Go is the solution of the equation 
where the 'velocity' components Vg are identically zero: (dt + DQk^)Go{k,t — t') = 6{t — t') 
which is: 

Go{k, t-t') = e{t - t')e-^ok\t-t') ^^Q^ 

Or, in real space, 



• Ambiguities due to multiplicative noise. Ito vs. Stratonovitch. 

In equation (^), we have omitted to specify the dependence on the variable t. There is 
actually an ambiguity in the product term WqVk-g- In the discrete Chicago model Q, the 
g-t's emitted from a given site are independent of the value of the weight on that site. In 
the continuum limit, this corresponds to choosing Wg{t) to be independent of Vk-g{t), or else 
that the f 's must be thought of as slightly posterior to the w's (i.e. the product is read as 
Wg{t — 0)vk-q{t + 0)). In this case, the average of equation (P) is trivial and coincides with 
the noiseless limit; hence G = Gq. This can be understood directly on the discrete model by 
noticing that the Green function G{i,j\0, 0) can be expressed as a sum over paths, all starting 
at site (0,0), and ending at site (i,i): 

G(i,j|0,0)= n ^±(^'0 (12) 

paths V (k,l)£P 



where the q±{k,l) are either q+{k,l) or q^{k,l), depending on the path. Since each bond 
q±{k, I) appears only once in the product, the averaging over q is trivial and leads to: 

G{i,j\0,0)= Yl 2"J^ = Go(i,j|0,0) (13) 
paths v 

(Note that this argument fails for the computation of the correlation function C, since paths 
can 'interfere'. We shall return later to this calculation.) 

The above choice corresponds to Ito's prescription in stochastic calculus. Another choice 
(i.e. Stratonovitch's prescription) is however possible, which corresponds to the proper con- 
tinuum time limit in the case where the correlation length if is very small, but not smaller 
than a (see figure ^). In this case, the w's and the u's cannot be taken to be independent. 
This is the choice that we shall make in the following. 




t-t' 



Figure 2: Correlation function of the noise along t axis. The results presented below would hold for 
an arbitrary symmetric, short range function. 



2.3 Calculation of the averaged response and correlation functions. 

Two approaches will be presented. The first one, based on Novikov's theorem, leads to exact 
differential equations for G and C, which can be fully solved. The second one is a mode- 
coupling approximation (mca), based on a resummation of perturbation theory. It happens 
that, for this particular model where the noise is gaussian and short range correlated in time, 
both approaches give the same results, because perturbation theory is trivial. In other cases, 
though, where exact solutions are no longer available, the MCA is in general very useful to 
obtain non perturbative results (see [p^]). 

We shall see that the effect of the noise is to widen the diffusion peak: Dq is renormalized 
by an additional term proportional to the variance of the noise v. 

• Novikov's theorem. Exact equations for G. 

Novikov's theorem provides the following identity, valid if the v are gaussian random 
variables: 

{wik,t)vik',t)) = l^dt'J dq(^^^j^yv{q,t')v{k',t)) (14) 
Such a term actually appears in equation (^, after transformation into an equation for G: 
{dt + Dok')G{k, t-t')= 6{t - t') - iky^^ I ^ {v{q, t)w{k - q, t)) (15) 



In the limit where ix = a ^ 0, the noise correlation is of the form: {v{q,t)v{q' ,t')) = 
27ra'^6{q + q')gx{q)gt{t - t'), with gt peaked in t = t' such that f{t')gt{t - t') ~ f{t)gt{t - t') 
for any function /. In all section ^ we take gx{q) = 1- From formally integrating equation (^) 
between t' and t, one can express the equal-time derivative 6w/6v as: 



5w{k, t) 



6v{k',t') 

and thus obtain: 



= -ikw{k- k',t) (16) 

t'=t-0 



{dt + Dok^)G{k, t-t') = 5{t - t') - 27ra^kG{k, t - t') /* dt'gt{t -t') f ^{k - q) (17) 

JO J 27r 

Using the shape of the function gt (see figure ^), the first integral is 1/2. The second one is a 
convolution integral, and its value is A.k/2'K + 0{k'^) (see equation (^)). The final differential 
equation for G is then, in the small-fc limit, a diffusion equation with a renormalized diffusion 
constant: 

Dr = Do + ^ (18) 

It is interesting to note that the model remains well defined in the limit where the 'bare' 
diffusion constant is zero, since a non zero diffusion constant is induced by the fluctuating 
velocity v. This would not be true if equation (^) was interpreted with the Ito convention, 
where the fluctuating velocity would not lead to any spreading of the average density. 

The most important conclusion is thus that, in the present scalar model, stresses prop- 
agates essentially vertically: taking £ ~ a, the response at depth H to a. small perturbation 
is confined within a distance oc ^JDrH from the vertical. Since Dr ~ ^^/a, \/DrH is much 
less than H in the limit where H ^> £^/a, i.e. when the height of the assembly of grains is 
much larger than the grain size. 

• Exact equations for C. 

Exact equations can also be derived for C, following very similar calculations. From equa- 
tion (P), one can deduce the corresponding one for w{k, t)w{—k, t). Upon averaging, Novikov's 
theorem has to be used on quantities such as {w{k,t)v{q,t)w{—k — q,t)), finally leading to: 

(dt + 2DRk')C{k, t) = a'k' J ^C{q, t) (19) 

One can formally integrate equation (19). It gives 

C{k, t) = Cik, 0)e-2^«'='* + f dt'e-^^^^^^*-*''^C{t') (20) 

where C{t') = J ^C{k,t'). Let us specify at this stage two specific initial conditions C{k,0) 
which can be of interest. We consider, for simplicity, a random packing of 'table tennis' balls 
with no mass (p = 0), but subject to a random overload of zero mean ({w{x,0)w{x' ,0)) = 
Aq6{x — x')) or to a constant overload {w{x,t = 0) = Bq). Therefore, C{k,0) = Aq in the 
first case, and C{k,0) = B^d^k) in the second one. Equation (^) is then solved in two steps: 
we first integrate it over k and find a closed equation for C, which can be solved in Laplace 
transform. We call E the conjugate variable of t. From C{E), we get C{t) and then finally 
compute C{x, t). 
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Figure 3: Correlation function for the case of a random overload. B has been rescaled by the factor 
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1/a-1/(87cD„t) 



Figure 4: Correlation function for the case of a uniform overload. B has been rescaled by the factor 

4[Do+£>K(7r-l)J • 



o Random overload: in the small E (large t) limit, we get C{E) ~ meaning 
C{t) = ao/\/t, with oq = Dp+DflfTr-i) v^s^Dfl " finally leads to the following expression for 
B{x,t) = C(0,t) - C(x,t) = l/2{[w{x,t) - w{{),t)f): 



B{x = 0,t) = 
B{xy^a,t) 



^0 



1 _ e sDnt _^ 
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2[Do + Dii{7T-l)]ya 8DRt 
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(211 



which is shown in figure |^. 

o Constant overload: in the same limit, we get C{E) ^ or C{t) = bo, where 60 



Dr 



2[Do+DRi-7T~l) 



Hence: 

B{x = 0,t) = 
B{x^O,t) ■■ 



^2 r2 
cr Bq 



i[Do + Djiin - 1)] 



1 1 - e 



(22) 



which has a form similar to that above: see figure |^. 

One could have performed the calculation with the Ito convention (corresponding to the 
Chicago model). The final results for the correlation function are actually very similar to 
those above. The main point is that the correlation is rather structureless. Equation (|2^ ) 
shows that the correlation function C{x > a,t) becomes zero for large times, a result that 
was used to establish equation (|2|). 

• Pertubation theory. 

The above method gives exact results, essentially because v{x,t) is short range correlated 
in time: Sw/6v is then only needed at coinciding times, where it is exactly known. This would 
not be true in general; furthermore Novikov's theorem requires v to be gaussian. It is thus 
interesting to show how a systematic perturbation scheme can be made to work by the use 
of diagrams to represent equation (|9|). The MCA (Mode Coupling Approximation) is then a 
particular resummation scheme of this set of diagrams, which was discussed in detail in [p^, 
which sometimes provide interesting non perturbative results. 

Equation (^) is multiplied on the left by the operator Go (see equation (p!o[)), and then 
reexpressed as follows: 



w{k,t) = Goik,t)(g) p{k,t) 



ikGo{k,t)^ I —■w{q,t)v{k - q,t) 



(23) 



(8) meaning a t-convolution product. This equation can be represented with diagrams as 
follows: as shown in figure |5|, we represent the source p by a cross, the 'bare' propagator Go 
by a plain line and the noise v by a dashed line. The first term of equation (p3|), which is the 
noiseless solution wq, is then obtained the juxtaposition of a plain line and a cross. The arrow 
flows against time (i.e. it is directed from t to t' < t). The juxtaposition of two objects means 
a t-convolution product. By definition w is represented by the juxtaposition of a bold line 
and a cross (this is consistent with the identification of a bold line with the full propagator 
G). The diagrammatic version of equation ([2^) is then: 



kk-q 

I 

k ' a 



p(k,t) = X 



Go(k,t) = 



k 
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-o- 



w{k, t) 



4< 



T.{k,t- t' 



k — q ^"'^ — k 



t 



t' 



Figure 5: definition of various diagrams. 



The 'vertex' stands for —ik J the two emerging wave vectors being q and k — q (node law). 
One can now iterate this equation. To second order, one obtains: 



^-K = > — X + ^ ' — ^ — X + > ' — > ' — > — X (25) 

The corresponding equation for G is obtained by taking the derivative S/6p, and averaging 
over the noise v. Since {v) = 0, the second diagram vanishes. We represent the noise correlator 
by a dashed line with a centered circle (see figure |5|), and obtain: 

/ 

* = ^ + ^ ' ^ ^ (26) 



or G = Go + Go S Go, where S is called the self-energy (see figure ^). Actually, one can resum 
exactly all the diagrams corresponding to GoSGq, GqSGoSGo to obtain the Dyson equation 
G = Go + GoSG. 

The MCA amounts to replacing the 'bare' propagator in the diagram for S by the full 
propagator G. (Note that the MCA is of course exact to second order in perturbation theory). 
We then obtain a self-consistent equation for G: 

1 1 ' \ 

5^MCA = Gq — G^jcA = ' * (27) 

Diagrams like the one drawn in figure ^ are now also included. The self-energy Smca can be 

^ \ 
/ \ 



Figure 6: Example of a diagram included in the MCA. 

easily computed, we get ^ 

^MCA{k,t- t') = -2na\ f ^qG^cAiq^t - t')gt{t - t') (28) 



In the special case where gt is peaked around t = t' , we can make the approximation G{q, t — 
t')gt{t-t') ~ G{q,0)gt{t-t') = gt{t-t') (since by definition G(g,0) = 1). We thus get, using 
equation (|), SMCA(^,i - t') = -a'^Ak'^gtit - t'). The expression for G^^^j^ is thus identical to 
the one obtained with the exact approach, as can be seen by comparing equation (0) and 

Note that one can also calculate the influence of a non zero kurtosis k of the noise v, 
which is its normalized fourth cumulant. In this case, four dashed lines (corresponding to v) 
can be merged, leading to a contribution to D, of the order of na^. 

Let us turn now to the calculation of the correlation function {w{k, t)w{k' , t)) = 2Tr6{k + 
k')C{k,t). The basic object which corresponds to the self-energy is now the 'renormalized 
source' spectrum S{k, t, t') defined as: C = G (8> S" (8> G. The quantity S is drawn as a filled 
square. 5*0 (empty square) is the correlation function source term which encodes the initial 
conditions (see below). The two first terms of the expansion are 



= ^ 




+ i r + ... 



(29) 

Here again, we transform the perturbative expansion into a closed self-consistent equation 
for S by replacing Go and So in (|2^) by G and S respectively. The final equation for C reads: 



^ ^ (30) 



or, written explicitly. 



C{k, t) = f dt' f dt"G{k, t - t')So{k, t',t")G{-k, t - t")+ 
Jo Jo 

a^k^ j dt' [ dt"G{k,t-t') [ ^G{q,t' ,t")gt{t' - t")G{-k,t - t") (31) 
JO Jo J 27r 

If we choose the source term to be an overload localised at t = 0, we get: Sq = {p{k, t')p{—k, t")) 
C{k,0)6{t')6{t"). 

Using the fact that gt is peaked around t' = t" , we again recover exactly the equation 
(|l9|) above, showing again that MCA is exact in this special case. 



2.4 Further results: the unaveraged response function 

The average Green function described above is thus a Gaussian of zero mean, and of width 
growing as \JD^t. However, for a given environment, the Green function is not Gaussian, 
presenting sample dependent peaks (see figure ^). Note however that, contrarily to what we 
shall find below for the tensorial case, the unaveraged Green function remains everywhere 
positive. Furthermore, the quantity [x](t), defined as the displacement of the centroid of the 
weight distribution beneath a point source in a given realisation: 




Figure 7: Averaged (bold line) and unaveraged (thin line) response function of the scalar model, 
obtained numerically by simulating the Chicago model. The average is performed over 5000 samples. 
One can notice how 'non self averaging' is the response function, i.e. how different it is for a given 
environment as compared to the average. Note also that the unaveraged Green function is everywhere 
positive. 



typically grows with t. More precisely, one can show that: 

{[x]{t))=0 but ([x]2(i)) oci^/^ (33) 

meaning that the 'center' of Green function wanders away from the origin in a subdiffusive 
fashion, as t^^'^. This behaviour has actually been obtained in an another context, that of 
a quantum particle interacting with a time dependent random environment. Physically, the 
Chicago model can indeed be seen as a collection of time dependent scatterers, converting 
ingoing waves into outgoing waves with a certain partition factor g+ = 1 — (see the 
discussion in [^]). In two dimensions (plus time), the wandering of the packet center [x]{t) 
is only logarithmic (and disappears in higher dimensions |l9| ). 

2.5 The scalar model with bias: Edwards' picture of arches 

Up to now, we have considered the mean value of v to be zero, which reflects the fact that 
there is no preferred direction for stress propagation. In some cases however, this may not be 
true. Consider for example a sandpile built from a point source: the history of the grains will 
certainly inprint a certain oriented 'texture' to the contact network, which can be modelled, 
within the present scalar model, as a non zero value of {v), the sign of which depends on 
which side of the pile is chosen. Let us call Vq the average value of v on the x > side of the 
pile, and — Vq on the other side. The differential equation describing propagation now reads, 
in the absence of disorder: 

dtw + dx [Vo sign{x)w] = p + Dod^xW (34) 
(An extra noise can be handled as above). For a constant density p = po, and for Dq = 0, 




0.5 

normalized radius 



Figure 8: On the main graph is plotted the solution of equation ( |34D for Vb/c = 0.4. The dashed line 
is for a diffusion constant Do ten times smaller than the the solid one. The bold line is for Dq = 0. 
Stress values are rescaled by the height of the pile t. The left inset graph shows that w{0, t) scales like 
at small Vb while the right inset shows that ^(Ojt) is constant for large t. Note that for very 
small values of Vb, the scaling becomes invalid for finite size reasons. 



the weight distribution is then the following: 



w{x^ t) 
w{x, t) 



Pox 
Vb 

po{ct-x) 
c-Vo 



for < X < Vot 
for Vot < X < ct 



(35) 



where c = l/tan(0) ((/> is the angle made by the slope of the pile with the horizontal x axis). 
For Dq 7^ 0, the above solution is smoothed (see figure Q). In any case, the local weight 
reaches a minimum around x = 0. Equation ( |3^ ) gives a precise mathematical content to 
Edwards' model of arching in sandpiles ||2^, as the physical mechanism leading to a 'dip' in 
the pressure distribution [^]. As discussed elsewhere [0, 12 1, this can be taken much further 
within a tensorial framework (see section |^ . 

The scaling of the stress at the center of the pile, w{0,t), can be understood simply in 
terms of random walks subject to a bias Vq. The region contributing to w{0,t) is then found 
to be of finite volume, independent of t, and of the order of Dq/Vq, as shown on the two top 
pictures of figure |8|. 

Equation (^4|) with noise can in fact be obtained naturally within an extended Chicago 
model, with an extra rule accounting for the fact that a grain can slide and lose contact with 
one of its two downward neighbours This generically leads to arching; in the sandpile 
geometry and for above a certain probability of (local) sliding, the effective 'velocity' Vq 
becomes non zero and the weight profile (35) is recovered |^]. However, this extra sliding rule 
implicitly refers to the existence of shear stresses, which are absent in the scalar model, but 
which are crucial to obtain symmetry breaking effects modelled by a non zero Vq. It is thus 
important to consider from the start the fact that stress has a tensorial, rather than scalar, 
nature. This is what we investigate in the following section. 



3 The Tensorial Model 



3.1 The wave equation 

It is useful to start with a simple 'toy' model for stress propagation, which is the analogue of 
the model presented in figure |l|. We now consider the case of three downward neighbours (see 
figure for a reason which will become clear below. Each grain transmits to its downward 




Figure 9: Three neighbour configuration. Each grain transmits two force components to its downward 
neighbours. A fraction p of the vertical component is transmitted through the middle leg, and a fraction 
{I— p)/2 through each of the external legs. 



neighbours not one, but two force components: one along the vertical axis t and one along 
x, which we call respectively Ft{i,j) and Fx{i,j)- (We will restrict attention, in the sequel, 
to two-dimensional piles, leaving extensions to three dimensions for further investigations.) 
For simplicity, we assume that the 'legs' emerging from a given grain can only transport the 
vector component of the force parallel to itself (but more general rules could be invented). 
Assuming that the transmission rules are locally symmetric, and that a fraction p < 1 of the 
vertical component travels through the middle leg, we find: 



Ft{i,3) 



^[F,(i-I,j-I)+F,(i + l,j-I)] 

+ _ p) tan V [Ft(i - 1, j - 1) - Ft{i + l,i - 1)] 

WQ+pFt{i,3 - 1) + \{^-p) [Ft{i - 1, j - 1) + Ft{i + 1,3- I)] 



1 



2tan^/' 



[F,(i-l,j-l)-F,(i + l,i-I)] 



(36) 



(37) 



where if) is the angle between grains, defined in figure Taking the continuum limit of the 
above equations leads to: 



dtFt + d^F^ 
dtFx + 



clFt 



P 




(38) 
(39) 



where Cg 



(1 — p) tan^ ■0. Eliminating (say) Fx between the above two equations leads to 
a wave equation for Ft, where the vertical coordinate t plays the role of time and cq is the 
equivalent of the 'speed of light'. In particular, the stress does not propagate vertically, as it 
does in the scalar model, but rather at a non zero angle ip such that cq = tan(^. Note that 
(/? 7^ in general (unless p = 0); the angle at which stress propagates has nothing to do with 
the underlying lattice structure and can in principle be arbitrary. We chose a three leg model 
to illustrate this particular point. 



The above derivation can be reformulated in terms of classical continuum mechanics as 
follows. Considering all stress tensor components aij, the equilibrium equation reads, 



Ot(7tt + OxCTxt 

XX 



p 





(40) 
(41) 



Identifying the local average of Ff with au and that of with crj^, we see that the above 
equations ( |38| , 39) are actually identical to ( ^ 41) provided atx = CTxt (which corresponds 
to the absence of local torque) and axx = c^c^tt- This relation between normal stresses was 



postulated in |10| as the simplest constitutive relation obeying the correct symmetries which 
enables one to lift the indeterminacy of equations ( ^ , ^); it can be seen as a local Janssen 
approximation jl^]. Cq should encode the relevant information of the local geometry of the 
packing, friction, shape of grains, etc., and should thus depend on the construction history 
of the grain assembly. For example, in the sandpile geometry Cg is related to the angle of 
friction (/> of the material by the relation Cq = 1/(1 + 2tan-^(/>) [^]. This approach can be 
generalized to take into account a local asymmetry in the packing texture (which one expects 
for example in the case of a sandpile constructed from a point source), by allowing Cq to 
depend on axt/(^tt [0, EH- dependence is linear, this is equivalent to a coordinate 

rotation in x, t 111 



3.2 A stochastic wave equation 

The starting point of the scalar model is thus essentially the diffusion equation, which one 
perturbs by adding a random convective term. As the above paragraph shows, a more natural 
starting point is the wave equation. The toy model presented above however suggests that, 
provided local conservation laws are obeyed (i.e. those arising from mechanical equilibrium). 



many local rules for force transmission are compatible with the contact conditions [15|. It is 
thus natural to encode the disorder of the packing, or model the indeterminacy of the contact 
conditions as a randomly varying 'speed of light' cq (reflecting the fact that, for example, the 
parameter p can vary from grain to grain) . Two recent numerical simulations ^ actually 



suggest that this should be a good first approximation. In figure 10, we show a scatter plot 
of Gxx versus au, measured as averages of the local forces over a small box centered around 
different points within a heap (from ref.|2^]). This plot clearly shows that a linear relation is 



indeed acceptable, leading in this case to Cq ~ 0.56 it 0.03 [23|. There are however significant 
fluctuations, reflecting some disorder in the packing, which are furthermore uncorrelated from 
point to point. The histogram of v defined as: 

(^xx = cl[l + v{x, t)]att {v{x, t)) = (42) 

is found to be roughly gaussian, of relative width a ^ 0.3. This corresponds to a locally 
varying angle of stress propagation, which varies around the mean angle 54° by an amount 
~ 10.8". 

Motivated by the simulations results, we now investigate a model (called 'random sym- 
metric model' in the following) with the inhomogeneous constitutive relation (^2[), which 
leads to the following stochastic wave equation for stress propagation: 

dttcrtt = dxx Co(l + v{x,t)) au (43) 

where the random noise v is assumed to be correlated as {v{x, t)v{x' ,t')) = a'^gx{x — x')gt{t — 
t'). The correlation lengths ix and it are again kept finite, and of the same order of mag- 
nitude. In Fourier transform, this relation can also be written {v{k,t)v{k' ,t')) = 2TTa'^6{k + 
k')gx(k)gt(t — t'). It turns out that the final shape of the averaged response function depends 







Figure 10: Relation between axx and au from a microscopic numerical simulation of grains forming 
a heap in two dimensions p3| . These data are compatible with a stochastic constitutive relation 
(Jxx — Cq[1 + v{x,t)]att, where u is a random noise. 



on the sign of gxi^)- In section |2| we implicitly made the choice gxik) = 1, which corresponds 
to: gx{x = 0) = 1/a and gx{x > 0) = 0. We will keep this choice for the following calculations, 
but note that another form for gx could lead to sign (^a; (A)) = — 1. 

In the following, au will be again denoted by w. After a Fourier transform along x-axis, 
we get, from equation (Hsl) 



{dtt + clk^)w = dtp 



dq 
2^ 



w{q,t)v{k - q,t) 



(44) 



Note that the 'source' term of this equation is now dtp rather than p itself. Therefore, if we call 
G the Green function (or propagator) of this equation G = {6w/5dtp); the response function 
R = {5w/5p) of our system is now actually the time derivative of G: R{k,t) = dtG{k,t). 

The noiseless propagator Gq is the solution of the ordinary wave equation {dtt+CQk'^)Go{k, t- 
t') = d{t — t') and can be easily calculated: 



GQ{k,t) 



1 



2icok 



Jcokt 



-icQkt 



e{t) 



which leads to the response function Rq 



1 



Mx. t) = i^[Hx- Cot) + 5{x + Cot)] e{t) 



(45) 



(46) 



This last equation sums up one of the major results of (see also |11, ^]): in two dimen- 
sions, stress propagates along two characteristic rays. (Note that the corresponding response 
function in three dimensions reads Ro{x,t) oc (cgt^ — x^)~^/^ for \x\ < coi and zero otherwise 



|10| ). A relevant question is now to ask how these rays survive in the presence of disorder. 
We will show that for weak disorder, the (5-peaks acquire a finite (diffusive) width, and that 



^In three dimensions a secondary closure is needed, for instance a^x = Oyy, y being the third coordinate. 



the 'speed of light' is renormaUzed to a lower value. Not surprisingly, the effect of disorder 
can be described by an 'optical index' n > 1. For a strong disorder, however, we find (within 
a gaussian approximation for the noise v) that the speed of light vanishes and then becomes 
imaginary. The 'propagative' nature of the stress transmission disappears and the system 
behaves more like an elastic body, in a sense clarified below. 



3.3 Calculation of the averaged response function 

One can again use Novikov's theorem in the present case if the noise is gaussian and short 
range correlated in time. However, the same results are again obtained within the diagram- 
matic approach explained in section |2|, which can be easily transposed to the present case, 
and is more general. The propagator G is a now represented as a line, the source dtp a cross 
and the vertex meaning — CgA;^ / Within the MCA, the self-consistent equation (analogous 



to equation (27) in the scalar case) is: 

rt 



{dtt + 4k^)H{k, t) = 5{t) + / dt'SMCA(fc, t')H{k, t - t') (47) 

JO 

where H is defined by G{k,t) = H{k,t)9{t), and the self-energy T^mca given as 

^MCA{k, t-t') = 27rcyk^ J ^q^9t{t - t')~g^{k - q)H{q, t - t') (48) 

Equation ( p7| ) can be solved using a standard Laplace transform along the t-axis {E is the 
Laplace variable). Using the fact that H{k,T) = r in the limit where r — > 0, we find, for 
small k, E (corresponding to scales L such that l^-, £t <C L): H~^{k, E) = E"^ + j3E + cj^k"^, 
where 



^2 4cT^A% (, m 



cm = ci-^-^^^^[l-^-^^+0{e) (49) 

„4 2 7,2 a3/'2 

m) = g ^* +0(fc^) (50) 

We notice here that in the limit It 0, the effect of the randomness completely disappears, 
as in the scalar model with the Ito convention. (Technically, this is due to the fact that 
G{k,t = 0) = in the present problem). In order to calculate the inverse Laplace transform, 
we need to know the roots of the equation H~^[k,E) = 0. This leads to several phases, 
depending on the strength of the disorder. 



• The weak disorder limit. 



For a weak disorder, cj^ik) is always positive. We can then define cji = CR{k = 0). As we 
will show now, cr is the shifted 'cone' angle along which stress propagates asymptotically. 
cji is a decreasing function of a, meaning that the peaks of the response function get closer 
together as the disorder increases 0. For a critical value|^ a = ac, cr vanishes, and becomes 
imaginary for stronger disorder. 

In the limit of large t, the propagator reads: 

G{k,t) = -\ sm[cRkt{l + a\k\)] e-^^'^e{t) (51) 

CRk 

*As a technical remark, let us note that if gt = gx, the problem is symmetric in the change x t, 
CQ{x,t) l/cQ{x,t). It thus looks as if the cone should both narrow or widen, depending on the arbitrary 
choice of x and t. There is however no contradiction with the above calculation since we assumed that v has 
zero mean, while 1/(1 + v) — 1 has a positive mean value, of order a^. 

^For £t = tj, = 1 and Cq = 0.6 (corresponding to 4> = 30°), one finds ac 2^ 0.57. 



where the following constants have been introduced |^: 



a 



7 



(3{k) _ a^A^e 



(52) 

2k^- 18 

From equation ([5l|), the response function R, in the limit of small k and large t, is given by: 

R{k, t) = cos [cRkt{l + a\k\)] e'^^^^eit) (54) 

or in real space, 

where the scaling variables measuring distances relative to the two peaks, are defined by 
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Figure 11: Response function for weak disorder (ct/ctc ~ 0.32). The two curves have been rescaled by 

1 /2 

the factor 2 [47r|7|t] ' . The main graph shows the general double-peaked shape of the response of the 
system when subjected to a peaked overload at a; = 0, t — 0. The inset gives details the right-hand 
peak as a function of the scaling variable Note the asymmetry (for .gx(A) > 0), compatible with 
the results found in ll^. Note also that the curve becomes negative around ^_ = 2. 



X±CRt 

and where 7 = 7 — icna and b = e^^^^^ . $ is the standard error function. Figure 11 shows 
R as given by expression (^). Interestingly, this propagator not only has a finite diffusive 
width oc but is also asymmetric around its maxima. Surprisingly, and in sharp contrast 
to the scalar case discussed above, the response function becomes negative in certain intervals 

®Note that the sign of a is dictated by the sign of 5a; (A). 



(although its integral is of course equal to one because of weight conservation). This means 
that pushing on a given point actually reduces the downward pressure on certain points. 
This can be interpreted as some kind of arching: increasing the shear stress does affect the 
propagation of the vertical stress, and may indeed lead to a reduction in its local value which is 
redistributed elsewhere. As we shall see in section |, the unaveraged response function indeed 
takes negative (and rather large) values. This is a very significant result since it suggests 
that granular materials may be susceptible to rearrangement under extremely weak external 
perturbations. Suppose indeed that as a result of the perturbation, a grain receives to a 
negative force larger than the preexisting vertical pressure. This grain will then move and a 
local rearrangement of contacts will occur, inducing a variation of co(x,t) as to reduce the 
cause of the instability. Thus, the stochastic wave equation implicitly demands rules similar to 
those introduced in Q to describe extreme sensitivity to external perturbations in silos. The 
present model, which is purely static, does not say what to do when a local rearrangement 
occurs, but certainly suggests that small perturbations will induce such rearrangements. 

It is interesting to note that this response function was numerically measured in ref. [ p^ ; 
its shape is compatible with the above expression; in particular, the two peaks were found to 
be asymmetric with a longer 'tail' extending inwards, as we obtain here. Note however that 
for 5x(A) < 0, the shape of the peaks is reversed: the small dips are located inside the peaks 
and the longer tail extends outwards. This is actually what we obtain numerically in section 



• Shear response function. 



Equation ( ^0| ) provides a straightforward way to calculate the shear response function Rg 
in terms of R. Indeed, one has: ikRs{k, t) = 6(t) — dtR{k, t). We thus get, in the limit of small 
k and large t, 

Rs{k, t) = -icR sin [cRkt(l + a|A:|)] e'"'''^^e{t) (57) 

This shear response function is very similar to R, except that it is, as expected, an odd 
function of x. 



• Effective large scale equations. 



It is interesting to know of which differential equations the response functions R and 
Rg, are solutions. These effective equations can be interpreted as a coarse-grained (hydro- 
dynamical) description of the propagation of a stress perturbation which takes into account 
the average effect of the local disorder. One problem however comes from the presence of 
the 'dispersion' term a\k\, which corresponds to a non-local operator in real space. We thus 
neglect this term in the following discussion, but one should keep in mind that the effective 
equation are actually non local. In any case, the main features of the response functions 
(peaks centered around x = iciji {cr < cq) with a diffusive width oc \/t) are not lost when 
setting a = (except for the fact that the response function can become negative which is 
related to a 7^ 0). Effective equations can then be written in the large t limit as: 

dt {Satt) = 6p- {6a^t) (58) 
dt {Suxt) = -CrOo: {6att) + 2jdx^ (Sa^t) (59) 

That disorder generates the diffusion terms 2jdxx {^<7xt) is rather intuitive and had been 
guessed in |1C]. This terms can be seen the first term of a gradient expansion of to the 



constitutive equations which have the correct symmetry, i.e.: 

{^(^xx) = Cr {6att) - jdx {Saxt) (60) 



{5a tx) = {S(Jxt) 



(61) 



where the second equation is imposed by the absence of local torque. 

We have thus shown that the introduction of a small disorder in the local direction of 
propagation does not change radically the nature of stress propagation on large length scales, 
although the peaks in the response function acquire a diffusive width. These peaks acquire a 
width of the order of ^J^H (where H is the height of the pile), and are thus well separated 
in the limit where ^ 7. As we shall see now, this is no longer true if the disorder becomes 
strong. 



• Critical disorder: The wave/ diffusion transition. 



When the disorder is so strong that cr just vanishes, the roots of H^^(k,E) = change 
nature, and so does the response function R. The two peaks of the previous expression for 
R merge together, while the width becomes anomalously large (oc t^/^). In the asymptotic, 
large t, regime we obtain: 



R{k,t) = e{t)cos \x\kf/'^t 



(62) 



where the new constant A is defined by A = coy/3/2A and 7 = Cq^j/S. The physical response 




Figure 12: Response function for a critical disorder: = 0. 



function R is plotted in figure 12, for different values of as a function of the scaling variable 

X 



At2/3 • 



(63) 



On the scale t^/^, the double peak structure of R is still visible. However, note that the term 
e — 'yk'^t cannot be neglected, even for large t; this means that the response function is never 



really a function of ^ only, as clear from figure 12. Note that the response function again 
becomes negative for some values of ^. 



Strong disorder: The pseudoelastic regime. 



For larger disorder still, one finds, within the MCA (which is exact for a gaussian, uncorre- 
lated noise), that the renormalized value of Cq, c^, becomes negative. Upon a rescaling of x as 
X = x/{icji), the effective equation on {6att) then becomes, on large length scales, Poisson's 
equation: 

V2 {6au) = dt {6p) (64) 

which means that the stress propagation becomes somewhat similar to that in an elastic 
body, where stresses obey an elliptic equation of similar type [^^. In particular, the cone 
structure of stress propagation, which is associated with the underlying, hyperbolic, wave 
equation finally disappears; the average response to a localised perturbation becomes a broad 
'bump' of width comparable to the height of the pile. It is thus rather interesting to see 
that, within MCA, there is a phase transition from a 'wavelike' mode to a 'diffusive' mode 
of stress propagation; the observation of the 'cone' thus requires that the packing is not too 
disordered. Certainly for relatively ordered packings the cone exists and has been observed 
experimentally and numerically [|l5|. One should however add some remarks: 

~ It is possible that the above transition is an artefact, due to the fact that v is taken to 
be gaussian, which strictly speaking is not allowed, since the local value of Cq should always 
be positive. One can show for some other problems of the same type that a similar transition 
is artificially induced by the gaussian approximation when it cannot really exist on physical 
grounds. In this respect, it is interesting to note that the first non gaussian correction tends 
to increase cr, for negative kurtosis as might be expected for a bounded v distribution. 

- It should be noted that the predicted effective constitutive relation between horizontal 
and vertical normal stresses has a negative sign if c| < 0. This means that increasing the 
vertical stress should reduce the horizontal stress, which is only possible is the grains move. 
Hence, the region where < is probably impossible to reach physically: the system will 
rearrange spontaneously as to reduce the disorder, and to make > 0. Note however that, as 
already discussed above, the disorder which results of such a rearrangement might be strongly 
correlated, and correspond to an arching effect, as in [^J. 

3.4 The correlation function 

Coming back to the weak disorder case, the major problem for the direct observation of the 
'light cone' is the fact that the perturbation representing the point source should be small 
(otherwise the packing structure would changes in an inhomogeneous way, thereby affecting 
the value of Cq in a non uniform way), but large enough for the response to be detected. A 
better possibility, as we show now, could be to measure the correlation function of the stress 
field. We again consider the stress correlation function in the case where the mass of each 
grain is small (p = 0) and a random or a constant overload is applied on the top of the 
'silo'. With the new convention for the bar (for G) and the cross (for dtp), the self-consistent 
diagrammatic equation ( |30|) is strictly valid in the tensorial case. When writing it into its 
usual mathematical form, the only difference with the scalar model is that now the weight 
source term is w{k,0)5'{t), leading to So{k,t',t") = C{k,0)5' {t')5' {t"). 

The calculation of the correlation function is very similar to the scalar case. In order 
to carry out the calculations to the end, we have neglected the dispersion term a\k\ in the 
expressions for G and R. The analogue of equation (|20|) is now, for weak disorder, 

G{k, t) = G{k, 0) cos^ [cRkt] e'^^^^'V 

a^^e f dt'sin^ [cRk{t-t')\ e-2^'='(*-*')C(t') (65) 



The function C{t') = J ^C{k,t') is of identical form to the scalar case; only the expressions 
for ao (for the random overload) and 60 (for the constant overload) are different: 
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(66) 
(67) 



Knowing C{t'), C{x,t) can be computed from equation (|6^). For the case of a constant 
overload, the shape of the correlation function is very close to the one showed in figure |^ for 
the scalar model. The case of the random overload however is much more interesting since the 
fact that information travels along a cone of angle cji appears clearly: the correlation function 
presents two peaks. The first one is of course at x = 0, while the second is at x = 2cRt, which 
simply means that the two points at the bottom of the information cone share the same 
information coming from the apex of this cone [|. If we forget the second term of the right 
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Figure 13: Correlation function for the case of a random overload. Note the presence of a peak 
centered in x = 2c fit which reflects the fact that information in the tensorial model is travelling along 
a cone of angle of cj^. In the case of a fluctuating density in the bulk of the pile, one should integrate 
( |68| ) with respect to t. The result is plotted in the inset: the correlation reaches rapidly a first plateau, 
and then increases again to a higher value around x — 2cBt. The relative difference of height between 
the two plateaus decreases as t^^l"^ . 



hand of equation (|65| ) which is negligible compared to the first one at large t, we can see that 
the second peak of the correlation function has a width oc \ft and a height oc This 



approximation is actually equivalent to saying that the (linear) effective equations (58, [59| ) 
are sufficient to calculate the correlation function for large times. Other source terms, such 
as a fiuctuating density in the bulk of the pile, can thus be easily accommodated by linear 

'^In the absence of disorder, the correlation function consists of two 5 peaks, one at 2; = and the other at 
X — 2cot of half the amplitude. 



superposition. We have thus plotted on figure ^ the quantity B{x,t) = C{0,t) — C{x,t), 
omitting the second term in the r.h.s. of equation (^). Analyticahy, we have 



B(xA) 
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2 + 2e 27 - 2e 



87t 



(68) 



This result is of importance since the shape of this correlation function clearly differs from the 
corresponding one in the scalar model. Measuring carefully the averaged correlation function 
of a granular system could then confirm (or disprove) the presence of a light ray-like propa- 
gation. In this respect, it is interesting to plot the correlation function for three dimensional 
packings as well. This correlation function only depends on the radial distance r between the 




r=2c„t 



Figure 14: Correlation function for three dimensional disordered packings with a random overload, 
neglecting again the second term in equation (|65|). Note that, as in two dimensions, the correlation 
function exhibits a peak around r = 2cRt. 



two points, as is plotted in figure 14. We note that, much as in two dimensions, the correlation 
decreases sharply on the scale of a few grains, but increases again for distances of the order of 
the height of the pile. Note that a stress correlation function was actually recently measured 
in and found to be featureless, but on very short scales x < 5a, as compared to the height 
of the pile H ~ 100a. We thus expect the features of the correlation function to show up on 
much larger scales ~ 2crH. 



4 Generalized wave equations 



It is tempting to generalize equations (|38|, |39|) and write the most general linear equations 
governing the propagation of the forces which are compatible with the (local) conservation 
rules. These equations were first written by de Gennes [^]: 



dtFt + d^[ri'{x,t)F^ + fi'{x,t)Ft] = p 
dtF^ + d^[r]{x,t)Ft + fi{x,t)F^] = 



(69) 
(70) 



Note that the terms /i, /x' break the symmetry x —x. This is allowed locally and does 
not show up on large scales if their average is zero. Another possibility (but without noise), 
considered in detail in [12|, is that fJ,{x,t) changes sign with x, i.e.: fJ,{x,t) = ji sign(x), which 



describes the fact that the texture of a sandpile depends on which 'side' of the pile one is 
looking at. Interestingly, equations (^, ^) still lead to wave-like propagation, but now the 
bisector of the 'light cone' makes a non zero angle with the vertical (when fi or /i' are non 



zero). In other words, equations (|69|, 70) describe a situation where not only the opening 
angle of the cone can vary in space, but also its average orientation. 

The same analytical techniques as above can be still be used. We shall only discuss some 
special cases 0: 

o Let us first set = ^u' = and consider the case where r/' is random, and r] fixed (and 
equal to Cg). Taking r]'{x,t) = t/q (1 + v{x,t)) with the noise v as above, one finds that the 
renormalized value of rj' is: 

r]R = r]o\l Q j (71) 

Now, on large length scales, one must recover the continuum equilibrium equations for the 



stress tensor, equations (40, 41). The condition of zero torque requires that the stress tensor 
is symmetric, and thus one must set 77^ = 1, which imposes a relation between t/q and the 
amplitude of the noise a. Note that beyond a certain value of a, this relation can no longer 
be satisfied with a real tjq. This again means that the packing is unstable mechanically and 
will rearrange so as to reduce the disorder. 

o Another interesting class of models, which one can call '/x models', is such that: rj = 
CQ,r]' = 1, but fJ-{x,t) = cov{x,t) and /x' = or vice-versa. These two cases yield identical 
results, namely, in the large t limit: 

Rik,t) = cos{cQkt)e-^^^^e{t) (72) 
Rs{k,t) = -icQs\ii{cokt)e-^^''^e{t) (73) 

where 7 = — . Note that in these cases, the response peaks acquire a finite diffusive width 
oc but the angle of the information cone is unaffected by the disorder (i.e. cq is not 
renormalized). 

o Finally, there are special 'symmetry' conditions where the equations can be decoupled 
and reduced to two 'scalar' models. We will refer to this case as the 'double scalar' model. 
This occurs when fi = fi' = coVi{x,t) and r]' = t]/cq = 1 -|- V2{x,t) where vi,V2 are two 
different sets of noise. Let us define a± = cqFi ± Fx, x± = x ^ cot and v± = vi ±V2, we then 
obtain: 



dta+ = cop - Co dx+ [v+a+] (74) 
dta- = cop- Co dx_[v-cr^] (75) 

showing that 0"+ and (T_ decouple, each propagating along two rays, of 'velocity' ±co, plus 
a small noise v± which, as in the scalar case, generates a nonzero diffusion constant. The 
response functions for au and (Txt are thus again made of two diffusive peaks of width oc \/t, 
centered in x = itcot. The interest of this double scalar limit is that one can deduce simply 
the probability distribution of the stresses from the Chicago model. This is developed below. 
Note also that by construction, this special form of disorder does not lead to negative vertical 
stresses. 



®To lowest order in perturbation theory, the case where disorder in present in the four terms rj,ri' , fj,, fi' 
simultaneously is very simply obtained by adding the corrections induced by each term taken individually. 



5 Stress distribution within the tensorial model 



A physically relevant question is to know how local stresses are distributed. We have seen 
above that within a scalar approach, an exponential-like distribution (possibly of the type 
exp — w^, with P > 1) is expected One can wonder whether this exponential distribution 

survives within a tensorial description, and what happens for very small stresses w — > 0. 
Unfortunately, the full distribution can only be computed analytically for the 'double scalar' 
model; but numerical results have also been obtained for the random symmetric model, and 
are described below. 

5.1 The 'double-scalar' limit 

In the 'double scalar' limit, the histogram of the stress distribution is obtained trivially by 
noting that since a+ = wi and (T_ = W2 travel along different paths, they are independent 
random variables. Taking cq to be unity for simplicity, one thus finds: 

PM = J dwi J dw2P*{wi)P*{w2)6iau - + ^^ ) (76) 

P{<7.t) = J dwi J dw2P*{wi)P*{w2)6{a,t - ~^^ ) (77) 

where P* is the distribution of weight pertaining to the scalar case, which, as mentionned 
above, depends on the specific form of the local disorder and on the discretisation procedure. 
In the strong disorder case which leads to equation (^ [in the case = 2], we thus find 
that P{(ytt) is still decaying exponentially (it is actually a T distribution of parameter 2A), 
although its variance is reduced by a factor 2. For N = 2, one simply gets 

PM = (78) 
PM = (kxt| + ^)e-2|--l (79) 

The preexponential factor is therefore noticeable different from the prediction of equation 
i). 



5.2 Numerical histograms for the random symmetric model and open 
problems 

The numerical analysis of equations (^, 41), with a stochastic constitutive relation axx = 
r][l + v{x, t)]att is actually not an easy task, and the final results depend rather sensitively on 
the chosen discretisation. For example, a naive discretisation of the random wave equation 
leads to a non zero diffusive width even in the absence of disorder, and thus makes it hard 
to measure the 'true' response function, which should, in the absence of disorder, consists of 
two 5 peaks. However, it should be noted that such diffusive term (or order a) are actually 
expected physically - they indeed appear when equations (^, ^) are expanded to second 
order in the lattice spacing. We shall come back to this point below. 

The method we chose is the following. Starting with points regularly arranged at t = 0, we 
construct the network of characteristics (in the mathematical sense). To each point {x,t) is 
associated a 'speed of light' co(x, t) = cov'[l + ^(^^j t)] which determines the directions of lines 
which propagate the component of the stress parallel to that line, away from the point {x, t). 
The point (x,t) is then generated by two 'parents' points {x',t') and {x",t") as indicated 
on figure |l5|-(a). It sometimes happens that the cone from {x',t') is so wide that it cannot 
intersects with the one from {x",t") [see figure |l5|-(b)]. We then impose x = x" and t = t" . 




Figure 15: The left-hand picture (a) shows the construction rule of the characteristics network: the 
'child' point {x,t) is located at the intersection of the cones from the two 'parent' points {x',t') and 
{x",t"). When the cones do not intersect (b), we choose {x,t) and {x",t") to be coincident. 



This actually can be viewed as a local kind of arching: the point (x" , t") not only supports its 
'parent' neighbours but also its 'same generation' neighbour {x',t'). This method has several 
advantages. Its physical interpretation is very clear: points are 'grains' and characteristics are 
'stress paths'. Figure 16 shows an example of the network of those paths. We can see how 
stress paths actually merge together and generate arching. Furthermore, there is strictly no 
diffusion in the absence of disorder, i.e. the Green function is exactly given by the sum of two 
(5-peaks. 

Although the noise v has been implicitly considered to be gaussian throughout this paper, 
for numerical simplicity we chose the following algorithm for the calculation of v{x,t). At 
each site {x,t), a random angle 6 is uniformly chosen between — and A^. A controls the 
amplitude of the noise. We then set co(x, t) = cq tan (7r/4 + 6{x, t)); v and 6 are then related 
by v{x,t) = ■ However, since the lattice itself is generated by the disorder, the 

precise correlation function Qx of the v's is not well controlled in this numerical scheme. This 
is rather important since we showed in section ^ that the structure of gx influences the shape 
of the response function (it determines whether the negative part of the response lies on 
the inward or outside edge of the main peak). In fact, the structure of the peaks we obtain 
numerically is reversed compared to that of our analytical calculation: see figure 17. 

The numerical histogram of the force distribution at the bottom of a 'silo' computed 
within this numerical scheme immediately reveals some problems. Since the lattice becomes 
more distorted as 'time' grows, the numerical histogram of vertical forces keeps broadening 
and never reaches a stationary shape. Furthermore, there is a nonzero probability of observ- 
ing negative weights which is, as we pointed out already, a structural property of the wave 
equation with randomness. Clearly, from a physical point of view, this is unacceptable and 
an additional rule should be added if the weight becomes locally negative. Some physically 
motivated rules could be invented (much as in [^), but we do not want to pursue this here, 
and leave this for future investigations. 

In the present paper, we restrict to the case of a nonzero 'bare' diffusion constant which, 
as argued above, should exist on a physical basis. Numerically, we have implemented this is 
two different ways. 

o The first one corresponds to letting the above scheme run until some height and then 
start afresh with a regular lattice, where the forces are obtained by averaging over the nearest 
neighbours belonging to the 'old' lattice. This averaging procedure is clearly equivalent to a 
diffusion term. In this case, the numerical histograms do reach a stationary limit. We note 
that: 
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Figure 16: Stress path network for a perfodic sifo of width 100 a. This picture has been computed 
with A = 0.2. We have chosen periodic lateral boundary conditions. 
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Figure 17: The main graph shows the response function calculated numerically on a silo of width 
100 a with A = 0.2. The thin line is a typical response for a given realisation of the disorder. Note 
that it takes negative values. The bold line has been averaged over 5000 realisation of the disorder. 
The inset compares the averaged response peak with the one computed analytically, with a negative 
a. Note the negative part, as predicted by the theoretical calculation. 
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Figure 18: These curves show the histograms of the vertical normal stress w, from which negative 
values have been removed. They all have been computed for the three-leg model, with periodic 'silos' 
of width 1000 a. The three bold (solid, long-dashed and dashed) lines are results from silos where 
the amplitude of the noise is maximum {pM — 1)- The height of those silos is as indicated in the 
legends. On the contrary, the thin line represent a 1000 a x 1000 a silo where the amplitude of the 
noise is pm = 0.5 where it is nearly gaussian. Much like within the scalar model, P{w) shows an 
exponential tail for large values of w when the disorder is maximum, while it is better fitted by a 
stretched exponential, logP{w) ^ —w^ with /3 > 1, for smaller values oipM- 



- The total probability of negative forces is reduced when to is smaller. 

- For tn ^ CL, the histogram is very nearly gaussian around the average force. 

- For larger to, the tail of the probability distribution for large forces is of the form 
exp—w^, where 2 > /? > 1 (as found in [^), where l3 is decreasing towards one as to 
increases, of for increasing disorder. For tf) = W a and A = 0.1, we found (3 ~ 1.6. The 
small force region has a much larger weight than found within the scalar model, although the 
presence of negative forces prevents us from being conclusive in this region. 

o The second scheme consists in simulating directly the three-leg model introduced above, 
with a random p chosen between and pm- These scheme is thus very close in spirit to the 
Chicago model. Again, the local forces are not everywhere positive, and thus the small force 
region cannot be reliably studied. Nevertheless, the large force region, however, behaves much 
in the same way as in the Chicago model. In particular, as shown in figure 18, the tail of the 
distribution decays as exp —w^, with /? ~ 1 when pM = 1, and with /3 > 1 when pM < 1- 

More work is needed to understand the physical implications of the presence of negative 
forces and any relation this may have to the static avalanche phenomenon However, the 
above results show that the tail of the force distribution is only exponential in a 'strong 
disorder' limit, where local 'arches' (i.e. one grain entirely bearing on a single downward 
neighbour) has a non zero probability. 



6 Summary- Conclusion 



We have investigated in great detail the role of a local disorder in the propagation of stresses 
in granular media, both within a 'scalar' approach, where only one component of the stress 
tensor is retained, and within the full tensorial approach, using a simple linear closure scheme 
(called 'bcc' in |11, |l2|), motivated partly by numerical simulations, which leads to a wave- like 
equation for stress propagation. The main effect of this local disorder is, besides introducing 
a diffusion like term in the effective, large scale equations, to renormalise the opening of the 
angle of the characteristic 'light cone' for propagation of stress. Within a 'Mode-Coupling' 
approximation (mca) scheme (exact for uncorrelated gaussian noise), one finds that this 
angle vanishes for a critical disorder, beyond which stress propagates in a fundamentally 
different way (this regime might however not be physically relevant). The most striking 
difference between the scalar and tensorial approach, is the fact that the response function 
becomes negative in the latter case, which is a source of instability of the packing to external 
perturbations. For moderate disorder, the response function takes negative values of order 
one near the point where the perturbation is applied, and decays with distance. Hence, we 
expect this instability to occur near the point where the perturbation is applied; at least 
near the upper surface of a pile under gravity, the effect occurs for a stress perturbation as 
small as the weight of one grain, since this is sufficient to make the total local vertical stress 
negative ! 

Another difference which could be amenable to experimental verification is the structure 
of the correlation function, which gives direct information on how the information travels in 
the medium. Because of the analogy between the scalar model and passive scalar convection 
in turbulence, it is furthermore possible that higher moments of the correlation functions 
might reveal, in some circumstances, an intermittent behaviour. Finally, the exponential fall 
off of the local stress distribution at high values, first found within the scalar model, also 
holds within a tensorial approach, but requires large disorder. 

Several open points remain for further studies. First of all, we have only considered two 
dimensional packings. The extension to three dimensions is rather straightforward - although 
the structure of the response functions becomes inherently more complex in this case (see 
|10|), the main features discussed here (i.e, diffusive spreading and narrowing of the cone) 
are still valid. 

Finally, we have not been able to determine analytically the histogram for local stresses 
within the random BCC model. The major unsolved problem is the presence of negative forces, 
which induce a mechanical instability, and imposes that an extra rule should be added to 
the stochastic wave equation to determine how stress propagates. As emphasized above, we 
believe that this is a direct consequence of the tensorial nature of the problem and can be 
interpreted as a signature of "fragility" of the contact network, which is generically unstable 
to very small perturbations ||3|, [27[| . 
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